# Start with removing all objects from possible previous syntax runs.
rm(list = ls(all=TRUE))


# Open packages used in this script file.
library(lme4)
library(foreign)
library(lattice)
library(readr)
library(dplyr)
library(haven)
library(ggplot2)
library(ggridges)
library(ggthemes)
library(gghighlight)
library(ggforce)
library(lmtest)
library(broom)
library(sjPlot)
library(broom.mixed)


mydata <- read_dta("C:/Users/tsipma/surfdrive/Wetenschappelijke artikelen/Emilien-David-Take/Data/Data analysis.dta")

str(mydata)
attach(mydata)

table(DMP)
table(REF)
table(trust_rep)
table(sat_dem)
table(income)
table(regioncode)
table(unemrate_2021)
table(gdp_2021)

#construct measure that combines support for DMP & REF
mydata <- mydata %>%
  mutate(DI = rowMeans(select(., DMP, REF), na.rm = TRUE))
table(mydata$DI)

correlation <- cor(mydata$REF, mydata$DMP, use = "complete.obs")
correlation


#descriptives
means_by_country_ref <- mydata %>%
  group_by(countrycode) %>%
  summarize(mean_REF = mean(REF, na.rm = TRUE),
            se_REF = sd(REF, na.rm = TRUE) / sqrt(n()), .groups = "drop")
print(means_by_country_ref)
means_by_region <- mydata %>%
  group_by(regioncode) %>%
  summarize(mean_unem = mean(unemrate_2021, na.rm = TRUE),
            mean_gdp = mean(gdp_2021, na.rm = TRUE),
          mean_EQI = mean(EQI_2021, na.rm = TRUE))
print(means_by_region)
correlation_reg1 <- cor(means_by_region$mean_unem, means_by_region$mean_gdp, use = "complete.obs")
correlation_reg1


#define directory to save
directory <-  "XXXXX"
theme_set(theme_classic())


####regression models 


#Model 1: ind+GDP (presented in Figure 2, appendix 6 Table A)
model1_di <- lmer (DI ~  
                   + gdp_2021 
                   + capitalregion
                   +    income  
                   + trust_rep + efficacy 
                   + interest + lrsp 
                   + urban + edu_group + factor(employment)  + male + age
                   + factor(countrycode) + (1|regioncode), 
                   data = mydata, verbose=TRUE)
model1_dmp <- lmer (DMP ~  
                     + gdp_2021 
                   + capitalregion
                   +    income  
                   + trust_rep + efficacy 
                   + interest + lrsp 
                   + urban + edu_group + factor(employment)  + male + age
                   + factor(countrycode) + (1|regioncode), 
                   data = mydata, verbose=TRUE)
model1_ref <- lmer (REF ~  
                     + gdp_2021 
                   + capitalregion
                   +    income  
                   + trust_rep + efficacy 
                   + interest + lrsp 
                   + urban + edu_group + factor(employment)  + male + age
                   + factor(countrycode) + (1|regioncode), 
                   data = mydata, verbose=TRUE)
summary(model1_di)
summary(model1_dmp)
summary(model1_ref)


#Model 2: only individual (appendix 8)
model2_di <- lmer (DI ~  
                   + capitalregion
                   +    income  
                   + trust_rep + efficacy 
                   + interest + lrsp 
                   + urban + edu_group + factor(employment)  + male + age
                   + factor(countrycode) + (1|regioncode), 
                   data = mydata, verbose=TRUE)
model2_dmp <- lmer (DMP ~  
                    + capitalregion
                    +    income  
                    + trust_rep + efficacy 
                    + interest + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + factor(countrycode) + (1|regioncode), 
                    data = mydata, verbose=TRUE)

model2_ref <- lmer (REF ~  
                    + capitalregion
                    +    income  
                    + trust_rep + efficacy 
                    + interest + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + factor(countrycode) + (1|regioncode), 
                    data = mydata, verbose=TRUE)
summary(model2_di)
summary(model2_dmp)
summary(model2_ref)

#Model 3: only GDP (appendix 8)
model3_di <- lmer (DI ~  
                     + gdp_2021 
                   + capitalregion
                   + interest + lrsp 
                   + urban + edu_group + factor(employment)  + male + age
                   + factor(countrycode) + (1|regioncode), 
                   data = mydata, verbose=TRUE)
model3_dmp <- lmer (DMP ~  
                      + gdp_2021 
                    + capitalregion
                    + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + factor(countrycode) + (1|regioncode), 
                    data = mydata, verbose=TRUE)
model3_ref <- lmer (REF ~  
                      + gdp_2021 
                    + capitalregion + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + factor(countrycode) + (1|regioncode), 
                    data = mydata, verbose=TRUE)
summary(model3_di)
summary(model3_dmp)
summary(model3_ref)

##interaction gdp*income (presented in Figure 3, appendix 6 Table B)
#12 income * gdp
model12_di <- lmer (DI ~  
                    + sub_income  
                    + trust_rep + efficacy 
                    + gdp_2021
                    + capitalregion
                    + interest + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + gdp_2021:sub_income
                    + factor(countrycode)
                    + (1+sub_income|regioncode), 
                    data = mydata, verbose=TRUE)
model12_dmp <- lmer (DMP ~  
                      + sub_income  
                    + trust_rep + efficacy 
                    + gdp_2021 
                    + capitalregion
                    + interest + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + gdp_2021:sub_income
                    + factor(countrycode)
                    + (1+sub_income|regioncode), 
                    data = mydata, verbose=TRUE)
model12_ref <- lmer (REF ~  
                      + sub_income  
                    + trust_rep + efficacy 
                    + gdp_2021 
                    + capitalregion
                    + interest + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + gdp_2021:sub_income
                    + factor(countrycode)
                    + (1+sub_income|regioncode), 
                    data = mydata, verbose=TRUE)

summary(model12_di)
#control = lmerControl(optimizer = "bobyqa") -> same results
summary(model12_dmp)
summary(model12_ref)

#22 trust_rep * gdp (presented in Figure 4, appendix 6 Table C)
model22_di <- lmer (DI ~  
                      + sub_income  
                    + trust_rep + efficacy 
                    + gdp_2021 
                    + capitalregion
                    + interest + lrsp
                    + urban + edu_group + factor(employment)  + male + age
                    + gdp_2021:trust_rep
                    + factor(countrycode)
                    + (1+trust_rep|regioncode), 
                    data = mydata, verbose=TRUE)
model22_dmp <- lmer (DMP ~  
                       + sub_income  
                     + trust_rep + efficacy 
                     + gdp_2021 
                     + capitalregion
                     + interest + lrsp 
                     + urban + edu_group + factor(employment)  + male + age
                     + gdp_2021:trust_rep
                     + factor(countrycode)
                     + (1+trust_rep|regioncode), 
                     data = mydata, verbose=TRUE)
model22_ref <- lmer (REF ~  
                       + sub_income  
                     + trust_rep + efficacy 
                     + gdp_2021 
                     + capitalregion
                     + interest + lrsp 
                     + urban + edu_group + factor(employment)  + male + age
                     + gdp_2021:trust_rep
                     + factor(countrycode)
                     + (1+trust_rep|regioncode), 
                     data = mydata, verbose=TRUE)

summary(model22_di)
summary(model22_dmp)
summary(model22_ref)

#42  + interest  * gdp (presented in Figure 5, appendix 6 Table D)
model42_di <- lmer (DI ~  
                      + interest
                         + sub_income  
                    + trust_rep + efficacy 
                    + gdp_2021 
                    + capitalregion
                    + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + gdp_2021:interest 
                    + factor(countrycode)
                    + (1+interest |regioncode), 
                    data = mydata, verbose=TRUE)
model42_dmp <- lmer (DMP ~  
                       + interest
                        + sub_income  
                     + trust_rep + efficacy 
                     + gdp_2021 
                     + capitalregion
                      + lrsp 
                     + urban + edu_group + factor(employment)  + male + age
                     + gdp_2021:interest
                     + factor(countrycode)
                     + (1+interest|regioncode), 
                     data = mydata, verbose=TRUE)
model42_ref <- lmer (REF ~  
                       + interest
                     + sub_income  
                     + trust_rep + efficacy 
                     + gdp_2021 
                     + capitalregion
                     + lrsp 
                     + urban + edu_group + factor(employment)  + male + age
                     + gdp_2021:interest
                     + factor(countrycode)
                     + (1+interest|regioncode), 
                     data = mydata, verbose=TRUE)

summary(model42_di)
summary(model42_dmp)
summary(model42_ref)

#robustness

#including gdp_change (Appendix 8)
model0_di_ch <- lmer (DI ~  
                     + gdp_2021  + gdp_change21_19 
                   + capitalregion
                   +    income  
                   + trust_rep + efficacy 
                   + interest + lrsp 
                   + urban + edu_group + factor(employment)  + male + age
                   + factor(countrycode) + (1|regioncode), 
                   data = mydata, verbose=TRUE)
summary(model0_di_ch)

model0_dmp_ch <- lmer (DMP ~  
                      + gdp_2021  + gdp_change21_19 
                    + capitalregion
                    +    income  
                    + trust_rep + efficacy 
                    + interest + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + factor(countrycode) + (1|regioncode), 
                    data = mydata, verbose=TRUE)
summary(model0_dmp_ch)

model0_ref_ch <- lmer (REF ~  
                      + gdp_2021  + gdp_change21_19 
                    + capitalregion
                    +    income  
                    + trust_rep + efficacy 
                    + interest + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + factor(countrycode) + (1|regioncode), 
                    data = mydata, verbose=TRUE)
summary(model0_ref_ch)


model0_di_ch11 <- lmer (DI ~  
                        + gdp_2021  + gdp_change21_11
                      + capitalregion
                      +    income  
                      + trust_rep + efficacy 
                      + interest + lrsp 
                      + urban + edu_group + factor(employment)  + male + age
                      + factor(countrycode) + (1|regioncode), 
                      data = mydata, verbose=TRUE)
summary(model0_di_ch11)
model0_di_ch16 <- lmer (DI ~  
                          + gdp_2021  + gdp_change21_16
                        + capitalregion
                        +    income  
                        + trust_rep + efficacy 
                        + interest + lrsp 
                        + urban + edu_group + factor(employment)  + male + age
                        + factor(countrycode) + (1|regioncode), 
                        data = mydata, verbose=TRUE)
summary(model0_di_ch16)

model0_dmp_ch16 <- lmer (DMP ~  
                          + gdp_2021  + gdp_change21_16
                        + capitalregion
                        +    income  
                        + trust_rep + efficacy 
                        + interest + lrsp 
                        + urban + edu_group + factor(employment)  + male + age
                        + factor(countrycode) + (1|regioncode), 
                        data = mydata, verbose=TRUE)
summary(model0_dmp_ch16)
model0_ref_ch16 <- lmer (REF ~  
                          + gdp_2021  + gdp_change21_16
                        + capitalregion
                        +    income  
                        + trust_rep + efficacy 
                        + interest + lrsp 
                        + urban + edu_group + factor(employment)  + male + age
                        + factor(countrycode) + (1|regioncode), 
                        data = mydata, verbose=TRUE)
summary(model0_ref_ch16)


###interactions with 5 year change
##interaction gdp*income
#12 income * gdp
model12_di_ch16 <- lmer (DI ~  
                    + sub_income  
                    + trust_rep + efficacy 
                    + gdp_2021  + gdp_change21_16
                    + capitalregion
                    + interest + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + gdp_change21_16:sub_income
                    + factor(countrycode)
                    + (1+sub_income|regioncode), 
                    data = mydata, verbose=TRUE)

#22 trust_rep * gdp
model22_di_ch16 <- lmer (DI ~  
                      + sub_income  
                    + trust_rep + efficacy 
                    + gdp_2021 + gdp_change21_16
                    + capitalregion
                    + interest + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + gdp_change21_16:trust_rep
                    + factor(countrycode)
                    + (1+trust_rep|regioncode), 
                    data = mydata, verbose=TRUE)

#42 interest  * gdp
model42_di_ch16 <- lmer (DI ~ + interest  
                           + sub_income  
                         + trust_rep + efficacy 
                         + gdp_2021 + gdp_change21_16
                         + capitalregion
                         + lrsp 
                         + urban + edu_group + factor(employment)  + male + age
                         + gdp_change21_16:interest 
                         + factor(countrycode)
                         + (1+interest|regioncode), 
                         data = mydata, verbose=TRUE)


summary(model12_di_ch16)
summary(model22_di_ch16)
summary(model42_di_ch16)


######################
##robustness for RR
##25 feb 2025

###participation instead of interest (Appendix 9)

#42  + participation instead of interest   * gdp (presented)
model42_par_di_direct <- lmer (DI ~  
                          + participation
                        + sub_income  
                        + trust_rep + efficacy 
                        + gdp_2021 
                        + capitalregion
                        + lrsp 
                        + urban + edu_group + factor(employment)  + male + age
                        + gdp_2021
                        + factor(countrycode)
                        + (1+participation|regioncode), 
                        data = mydata, verbose=TRUE)
summary(model42_par_di_direct)
model42_par_di <- lmer (DI ~  
                      + participation
                    + sub_income  
                    + trust_rep + efficacy 
                    + gdp_2021 
                    + capitalregion
                    + lrsp 
                    + urban + edu_group + factor(employment)  + male + age
                    + gdp_2021:participation
                    + factor(countrycode)
                    + (1+participation|regioncode), 
                    data = mydata, verbose=TRUE)
summary(model42_par_di)
model42_par_dmp <- lmer (DMP ~  
                       + participation
                     + sub_income  
                     + trust_rep + efficacy 
                     + gdp_2021 
                     + capitalregion
                     + lrsp 
                     + urban + edu_group + factor(employment)  + male + age
                     + gdp_2021:participation
                     + factor(countrycode)
                     + (1+participation|regioncode), 
                     data = mydata, verbose=TRUE)
model42_par_ref <- lmer (REF ~  
                       + participation
                     + sub_income  
                     + trust_rep + efficacy 
                     + gdp_2021 
                     + capitalregion
                     + lrsp 
                     + urban + edu_group + factor(employment)  + male + age
                     + gdp_2021:participation
                     + factor(countrycode)
                     + (1+participation|regioncode), 
                     data = mydata, verbose=TRUE)

summary(model42_par_di_direct)
summary(model42_par_di)
summary(model42_par_dmp)
summary(model42_par_ref)



#####FIGURES + Tables###############################################################


#tables
#to present
tab_model(model1_di, model1_dmp, model1_ref, show.se = TRUE, show.ci = FALSE, p.style = "stars")

tab_model(model12_di, model12_dmp, model12_ref, show.se = TRUE, show.ci = FALSE, p.style = "stars")
tab_model(model22_di, model22_dmp, model22_ref, show.se = TRUE, show.ci = FALSE, p.style = "stars")
tab_model(model42_di, model42_dmp, model42_ref, show.se = TRUE, show.ci = FALSE, p.style = "stars")

tab_model(model1_di, model2_di, model3_di, show.se = TRUE, show.ci = FALSE, p.style = "stars")

#change (appendix 8)
tab_model(model0_di_ch, model0_di_ch16, model0_di_ch11,  show.se = TRUE, show.ci = FALSE, p.style = "stars")
summary(model0_di_ch)

#political participation (appendix 9)
tab_model(model42_par_di, model42_par_dmp, model42_par_ref, show.se = TRUE, show.ci = FALSE, p.style = "stars")



#Figures for interaction
#model12 - income (figure 3)
fig1_di <- plot_model(model12_di, type="int", mdrt.values="minmax", 
                      title = "", axis.title = "Support for DIs",
                      legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Income security") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig1_di)
fig1_dmp <- plot_model(model12_dmp, type="int", mdrt.values="minmax", 
                      title = "", axis.title = "Support for DMP",
                      legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Income security") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig1_dmp)
fig1_ref <- plot_model(model12_ref, type="int", mdrt.values="minmax", 
                      title = "", axis.title = "Support for referendum",
                      legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Income security") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig1_ref)

fig1_class_di <- plot_model(model12_class_di, type="int", mdrt.values="minmax", 
                      title = "", axis.title = "Support for DIs",
                      legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Occupational class") + 
  ylim(c(1,5)) +
  scale_x_continuous(breaks = c(1,2), labels = c("Lower classes", "Higher classes")) +  
  theme(legend.position = "bottom")
print(fig1_class_di)
fig1_workingclass_di <- plot_model(model12_workingclass_di, type="int", mdrt.values="minmax", 
                            title = "", axis.title = "Support for DIs",
                            legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Production workers vs other occupational classes") + 
  ylim(c(1,5)) +
  scale_x_continuous(breaks = c(0,1), labels = c("Production workers", "Other classes")) +  
  theme(legend.position = "bottom")
print(fig1_workingclass_di)


#model22 - trust  (figure 4)
fig2_di <- plot_model(model22_di, type="int", mdrt.values="minmax", 
                      title = "", axis.title = "Support for DIs",
                      legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Trust in representative institutions") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig2_di)
fig2_dmp <- plot_model(model22_dmp, type="int", mdrt.values="minmax", 
                       title = "", axis.title = "Support for DMP",
                       legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Trust in representative institutions") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig2_dmp)
fig2_ref <- plot_model(model22_ref, type="int", mdrt.values="minmax", 
                       title = "", axis.title = "Support for referendum",
                       legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Trust in representative institutions") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig2_ref)


#model42 - interest  (figure 5)
fig4_di <- plot_model(model42_di, type="int", mdrt.values="minmax", 
                      title = "", axis.title = "Support for DIs",
                      legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Interest in politics") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig4_di)
fig4_dmp <- plot_model(model42_dmp, type="int", mdrt.values="minmax", 
                       title = "", axis.title = "Support for DMP",
                       legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Interest in politics") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig4_dmp)
fig4_ref <- plot_model(model42_ref, type="int", mdrt.values="minmax", 
                       title = "", axis.title = "Support for referendum",
                       legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Interest in politics") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig4_ref)

#model42 - participation instead of interest (appendix 9)
fig4_par_di <- plot_model(model42_par_di, type="int", mdrt.values="minmax", 
                      title = "", axis.title = "Support for DIs",
                      legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Political participation") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig4_par_di)
fig4_par_dmp <- plot_model(model42_par_dmp, type="int", mdrt.values="minmax", 
                       title = "", axis.title = "Support for DMP",
                       legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Political participation") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig4_par_dmp)
fig4_par_ref <- plot_model(model42_par_ref, type="int", mdrt.values="minmax", 
                       title = "", axis.title = "Support for referendum",
                       legend.title = "GDP per capita (in 10,000)",   colors = c("#003d5a","#ce132c")) +
  xlab("Political participation") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig4_par_ref)



##save figures

ggsave(file.path(directory, "Fig1.png"), fig1_di, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig1 dmp.png"), fig1_dmp, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig1 ref.png"), fig1_ref, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig1 workingclass.png"), fig1_workingclass_di, width = 4,
       height = 4, dpi = 300)

ggsave(file.path(directory, "Fig2.png"), fig2_di, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig2 dmp.png"), fig2_dmp, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig2 ref.png"), fig2_ref, width = 4,
       height = 4, dpi = 300)

ggsave(file.path(directory, "Fig3.png"), fig3_di, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig3 dmp.png"), fig3_dmp, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig3 ref.png"), fig3_ref, width = 4,
       height = 4, dpi = 300)

ggsave(file.path(directory, "Fig4.png"), fig4_di, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig4 dmp.png"), fig4_dmp, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig4 ref.png"), fig4_ref, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig4 participation.png"), fig4_par_di, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig4 participation dmp.png"), fig4_par_dmp, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig4 participation ref.png"), fig4_par_ref, width = 4,
       height = 4, dpi = 300)

ggsave(file.path(directory, "Fig5.png"), fig5_di, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig5 dmp.png"), fig5_dmp, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig5 ref.png"), fig5_ref, width = 4,
       height = 4, dpi = 300)

#coefplot for main text (figure 2)
#based on model 1


#coef model
#use model 1
#update 8-8-2024: efficacy deleted from model
model1_di_tidy <- tidy(model1_di, conf.int = TRUE)
model1_dmp_tidy <- tidy(model1_dmp, conf.int = TRUE)
model1_ref_tidy <- tidy(model1_ref, conf.int = TRUE)


model1_di_tidy <- model1_di_tidy  %>% mutate(model = c("Model 1"))
model1_dmp_tidy <- model1_dmp_tidy  %>% mutate(model = c("Model 2"))
model1_ref_tidy <- model1_ref_tidy  %>% mutate(model = c("Model 3"))

model1_di_tidy <- model1_di_tidy  %>% mutate(dv = c("1 Democratic innovations"))
model1_dmp_tidy <- model1_dmp_tidy  %>% mutate(dv = c("2 Deliberative mini-publics"))
model1_ref_tidy <- model1_ref_tidy  %>% mutate(dv = c("3 Referenda"))

model_total <- rbind(model1_di_tidy, model1_dmp_tidy, model1_ref_tidy)

fig_model_total <- subset(model_total, term == "gdp_2021" 
                          | term ==  "income" | term ==  "trust_rep" | term ==  "interest")

fig_model_total$order <- ifelse(fig_model_total$term == "gdp_2021", 1, ifelse(fig_model_total$term == "income", 2,
                                                                                     ifelse(fig_model_total$term == "trust_rep", 3,
                                                                                            ifelse(fig_model_total$term == "interest", 4, 99))))
#define limit
#limit
ylimit <- c(-0.3, 0.3)

#significant
#model25  <- model25  %>% 
#  mutate(significant = as.factor(ifelse(p.value < 0.05, 1, 0))) %>% 
#  filter(term != "(Intercept)")

fig_total <- ggplot(fig_model_total, aes(x = reorder(term, -order), y = estimate)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(size=2) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1), linewidth=0.5) +
  scale_x_discrete(labels = c("Individual political engagement (interest in politics)",  
                              "Individual political disaffection (trust in representative institutions)",
                              "Individual economic deprivation (income security)",
                               "Regional economic productivity (GDP/capita, in 10,000)")) +
  ylim(ylimit) +
  labs(x = "", y = "") +
  coord_flip() +
  facet_wrap(~ dv, ) +
  theme_classic() +
  theme(text = element_text(size = 18)) +
  theme(axis.text.x = element_text(size = 12),
        panel.spacing.x = unit(0.2, "cm"))

print(fig_total)
#ggsave(file.path(directory, "Fig_coef old.png"), fig_total, width = 12,
#       height = 4, dpi = 300)

ggsave(file.path(directory, "Fig_coef 5.png"), fig_total, width = 14,
       height = 5, dpi = 300)
ggsave(file.path(directory, "Fig_coef.png"), fig_total, width = 14,
       height = 4, dpi = 300)

#old labels
fig_total_oldlabel <- ggplot(fig_model_total, aes(x = reorder(term, -order), y = estimate)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(size=2) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1), linewidth=0.5) +
  scale_x_discrete(labels = c("Interest in politics",  
                              "Trust in representative institutions",
                              "Income security",
                              "GDP per capita (in 10,000)")) +
  ylim(ylimit) +
  labs(x = "", y = "") +
  coord_flip() +
  facet_wrap(~ dv, ) +
  theme_classic() +
  theme(text = element_text(size = 18)) +
  theme(axis.text.x = element_text(size = 12),
        panel.spacing.x = unit(0.2, "cm"))

print(fig_total_oldlabel)
ggsave(file.path(directory, "Fig_coef old.png"), fig_total_oldlabel, width = 12,
       height = 4, dpi = 300)

#coef model with efficacy
model1_di_tidy <- tidy(model1_di, conf.int = TRUE)
model1_dmp_tidy <- tidy(model1_dmp, conf.int = TRUE)
model1_ref_tidy <- tidy(model1_ref, conf.int = TRUE)


model1_di_tidy <- model1_di_tidy  %>% mutate(model = c("Model 1"))
model1_dmp_tidy <- model1_dmp_tidy  %>% mutate(model = c("Model 2"))
model1_ref_tidy <- model1_ref_tidy  %>% mutate(model = c("Model 3"))

model1_di_tidy <- model1_di_tidy  %>% mutate(dv = c("1 Democratic innovations"))
model1_dmp_tidy <- model1_dmp_tidy  %>% mutate(dv = c("2 Deliberative mini-publics"))
model1_ref_tidy <- model1_ref_tidy  %>% mutate(dv = c("3 Referenda"))

model_total <- rbind(model1_di_tidy, model1_dmp_tidy, model1_ref_tidy)

fig_model_total <- subset(model_total, term == "gdp_2021" | term ==  "capitalregion" 
                          | term ==  "income" | term ==  "trust_rep"
                          | term ==  "efficacy" | term ==  "interest")

fig_model_total$order <- ifelse(fig_model_total$term == "gdp_2021", 1, ifelse(fig_model_total$term == "capitalregion", 2, 
                                                                                   ifelse(fig_model_total$term == "income", 3,
                                                                                          ifelse(fig_model_total$term == "trust_rep", 4,
                                                                                                 ifelse(fig_model_total$term == "efficacy", 5,     
                                                                                                        ifelse(fig_model_total$term == "interest", 6, 99))))))
#define limit
#limit
ylimit <- c(-0.3, 0.3)

#significant
#model25  <- model25  %>% 
#  mutate(significant = as.factor(ifelse(p.value < 0.05, 1, 0))) %>% 
#  filter(term != "(Intercept)")

fig_total <- ggplot(fig_model_total, aes(x = reorder(term, -order), y = estimate)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(size=2) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0.1), linewidth=0.5) +
  scale_x_discrete(labels = c("Political interest", "Internal political efficacy", 
                              "Political trust",
                              "Subjective income",
                              "Capital region",
                              "GDP per capita (in 10,000)")) +
  ylim(ylimit) +
  labs(x = "", y = "") +
  coord_flip() +
  facet_wrap(~ dv, ) +
  theme_classic() +
  theme(text = element_text(size = 18)) +
  theme(axis.text.x = element_text(size = 12),
        panel.spacing.x = unit(0.2, "cm"))

print(fig_total)
ggsave(file.path(directory, "Fig_coef_efficacy.png"), fig_total, width = 12,
       height = 4, dpi = 300)


##robustness: change interactions (see appendix 8)

#model12 - income *change
fig1_di_ch16 <- plot_model(model12_di_ch16, type="int", mdrt.values="minmax", 
                      title = "", axis.title = "Support for DIs",
                      legend.title = "5-year change in GDP per capita",   colors = c("#003d5a","#ce132c")) +
  xlab("Income security") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig1_di_ch16)
#model22 - trust  *change
fig2_di_ch16 <- plot_model(model22_di_ch16, type="int", mdrt.values="minmax", 
                      title = "", axis.title = "Support for DIs",
                      legend.title = "5-year change in GDP per capita",   colors = c("#003d5a","#ce132c")) +
  xlab("Trust in representative institutions") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig2_di_ch16)
#model32 - efficacy  *change
fig3_di_ch16 <- plot_model(model32_di_ch16, type="int", mdrt.values="minmax", 
                      title = "", axis.title = "Support for DIs",
                      legend.title = "5-year change in GDP per capita",   colors = c("#003d5a","#ce132c")) +
  xlab("Internal political efficacy") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig3_di_ch16)
#model42 - interest  *change
fig4_di_ch16 <- plot_model(model42_di_ch16, type="int", mdrt.values="minmax", 
                      title = "", axis.title = "Support for DIs",
                      legend.title = "5-year change in GDP per capita",   colors = c("#003d5a","#ce132c")) +
  xlab("Interest in politics") + 
  ylim(c(1,5)) +
  theme(legend.position = "bottom")
print(fig4_di_ch16)

ggsave(file.path(directory, "Fig1 change.png"), fig1_di_ch16, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig2 change.png"), fig2_di_ch16, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig3 change.png"), fig3_di_ch16, width = 4,
       height = 4, dpi = 300)
ggsave(file.path(directory, "Fig4 change.png"), fig4_di_ch16, width = 4,
       height = 4, dpi = 300)

